# Replication Material for "Social Capital, Institutional Rules 
#     and Constitutional Amendment Rates"
# William Blake, Joseph Cozza, Dave Armstrong and Amanda Friesen
#
# before running, execute setup.r first
#
# This file produces tables Tables A6 and A7 and Figures 2, A6 and A7. 


cnl_dat <- import("crossnational_longitudinal_replication.dta")
cnl_ctrl_dat <- import("crossnational_longitudinal_ctrls_replication.dta")

m1b <- brm(amendment ~ vp3_b + words + cspi_w + tv + vp3_w:tv + 
             (1 | cowcode) + (1|spell), 
           data=filter(cnl_dat, lag_amend == 0), 
           family=bernoulli(),
           control = list(adapt_delta=.9), 
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars  = save_pars(all=TRUE))

m2b <- brm(amendment ~ vp3_b + words + cspi_w  +  vp3_w + 
             (1 | cowcode) + (1|spell), 
           data=filter(cnl_dat, lag_amend == 0 & tv == 1), 
           family=bernoulli(),
           control = list(adapt_delta=.99), 
           cores = 4, 
           backend = "cmdstanr", 
           refresh=0, 
           silent=2, 
           seed=207, 
           save_pars  = save_pars(all=TRUE))


m1bc <- brm(amendment ~ vp3_b + words + cspi_w + tv + vp3_w:tv + 
              exconst_w + gdppc_w + gdp_growth_w + war_w + EFindex_w + 
              (1 | cowcode) + (1|spell), 
            data=filter(cnl_ctrl_dat, lag_amend == 0), 
            family=bernoulli(),
            control = list(adapt_delta=.9), 
            cores = 4, 
            backend = "cmdstanr", 
            refresh=0, 
            silent=2, 
            seed=207, 
            save_pars  = save_pars(all=TRUE))

m2bc <- brm(amendment ~ vp3_b + words + cspi_w  +  vp3_w + 
              exconst_w + gdppc_w + gdp_growth_w +  war_w + EFindex_w + 
              (1 | cowcode) + (1|spell), 
            data=filter(cnl_ctrl_dat, lag_amend == 0 & tv == 1), 
            family=bernoulli(),
            control = list(adapt_delta=.99), 
            cores = 4, 
            backend = "cmdstanr", 
            refresh=0, 
            silent=2, 
            seed=207, 
            save_pars  = save_pars(all=TRUE))


## Table A6
s1f <- summary(m1b)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,4,5,2,6,3), 
         p = unname(c(sigf(m1b)))) %>% 
  arrange(obs) %>% 
  select(-obs)

s1r_c <- summary(m1b)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m1b, "cowcode"))))

s1r_s <- summary(m1b)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m1b, "spell"))))

s1 <- bind_rows(s1f, s1r_c, s1r_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") 

s2f <- summary(m2b)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(p = unname(c(sigf(m2b)))) 

s2r_c <- summary(m2b)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m2b, "cowcode"))))

s2r_s <- summary(m2b)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m2b, "spell"))))

s2 <- bind_rows(s2f, s2r_c, s2r_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") 

s1$term <- ifelse(s1$term == "tv:vp3_w", "vp3_w", s1$term)


s1 <- left_join(s1, s2 %>% rename(vals2 = vals)) %>% 
  select(-param)

s1$term[seq(2,16,by=2)] <- ""

s1 <-s1 %>% 
  mutate(term = gsub("vp3", "Const. Rigidity", term), 
         term = gsub("tv", "Time Varying", term), 
         term = gsub("cspi", "Soc. Capital", term), 
         term = gsub("judrev", "Jud. Review", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("sigma_", "Sigma: ", term), 
         term = gsub("_w", " (within)", term), 
         term = gsub("_b", " (between)", term))

s1$vals2 <- ifelse(is.na(s1$vals2), "", s1$vals2)
s1 <- setNames(s1, c("Term", "Model 1", "Model 2"))

x0 <- cnl_dat %>% filter(lag_amend == 0)
x1 <- cnl_dat %>% filter(lag_amend == 0 & tv == 1)
dim(x0)
dim(x1)
length(unique(x0$cowcode))
length(unique(x1$cowcode))
aux <- data.frame(Term = c("N", "N Countries"), 
                  `Model 1` = c("2002", "80"), 
                  `Model 2` = c("569", "16"))
names(aux) <- names(s1)
s1 <- bind_rows(s1, aux)

s1fc <- summary(m1bc)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,4,5,2,6,7,8,9,10,11,3), 
         p = unname(c(sigf(m1bc)))) %>% 
  arrange(obs) %>% 
  select(-obs)

s1rc_c <- summary(m1bc)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m1bc, "cowcode"))))

s1rc_s <- summary(m1bc)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m1bc, "spell"))))

s1c <- bind_rows(s1fc, s1rc_c, s1rc_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="valsc") 

s2fc <- summary(m2bc)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(p = unname(c(sigf(m2bc)))) 

s2rc_c <- summary(m2bc)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m2bc, "cowcode"))))

s2rc_s <- summary(m2bc)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m2bc, "spell"))))

s2c <- bind_rows(s2fc, s2rc_c, s2rc_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="valsc") 

s1c$term <- ifelse(s1c$term == "tv:vp3_w", "vp3_w", s1c$term)


s1c <- left_join(s1c, s2c %>% rename(vals2c = valsc)) %>% 
  select(-param)

s1c$term[seq(2,26,by=2)] <- ""

s1c <-s1c %>% 
  mutate(term = gsub("vp3", "Const. Rigidity", term), 
         term = gsub("tv", "Time Varying", term), 
         term = gsub("cspi", "Soc. Capital", term), 
         term = gsub("judrev", "Jud. Review", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("EFindex", "Ethnic Frac.", term), 
         term = gsub("exconst", "Exec. Const.", term), 
         term = gsub("war", "Conflict", term), 
         term = gsub("gdppc", "log(GDP/capita)", term), 
         term = gsub("gdp_growth", "GDP Growth", term), 
         term = gsub("sigma_", "Sigma: ", term), 
         term = gsub("_w", " (within)", term), 
         term = gsub("_b", " (between)", term))

s1c$vals2c <- ifelse(is.na(s1c$vals2c), "", s1c$vals2c)
s1c <- setNames(s1c, c("Term", "Model 3", "Model 4"))

z0 <- cnl_ctrl_dat %>% filter(lag_amend == 0)
z1 <- cnl_ctrl_dat %>% filter(lag_amend == 0 & tv == 1)
nrow(z0)
nrow(z1)
length(unique(z0$cowcode))
length(unique(z1$cowcode))

auxc <- data.frame(Term = c("N", "N Countries"), 
                   `Model 1` = c("1529", "78"), 
                   `Model 2` = c("397", "16"))
names(auxc) <- names(s1c)
s1c <- bind_rows(s1c, auxc)

blnk_row <- c("", "", "")
s1new <- rbind(s1[1:12, ], blnk_row, blnk_row, blnk_row, blnk_row, blnk_row, 
               blnk_row, blnk_row, blnk_row, blnk_row, blnk_row, s1[13:18, ])
s1_all <- bind_cols(s1c, s1new %>% select(-1)) %>% select(1,4,5,2,3)


l1 <- loo(m1b, moment_match=TRUE)
l2 <- loo(m2b, moment_match=TRUE)
l3 <- loo(m1bc, moment_match=TRUE)
l4 <- loo(m2bc, moment_match=TRUE)

l_ests <- rbind(l1$estimates[3,], 
                l2$estimates[3,],
                l3$estimates[3,],
                l4$estimates[3,])

loos <- c("LOO IC (SE)", apply(l_ests, 1, function(x)sprintf("%.1f (%.1f)", x[1], x[2])))
names(loos) <- names(s1_all)
s1_all <- rbind(s1_all, loos)

prefun <- function(obj){
  pcp <- mean(round(fitted(obj)[,1]) == obj$data$amendment)
  pmc <- mean(obj$data$amendment)
  pmc <- ifelse(pmc < .5, 1-pmc, pmc)
  (pcp-pmc)/(1-pmc)
}

r <- c("PRE", sprintf("%.2f", c(prefun(m1b), prefun(m2b), prefun(m1bc), prefun(m2bc))))
names(r) <- names(s1_all)
s1_all <- rbind(s1_all, r)
s1_all


## Figure A3

variable <- "cspi_w"
D1 <- D2 <- cnl_dat %>% 
  mutate(tv = as.numeric(tv))
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- -.5*s
D2[[variable]] <- .5*s

fit1 <- fitted(m1b, 
               newdata=D1,
               #re_formula=NA, 
               allow_new_levels=TRUE, 
               #                   sample_new_levels = "old_levels", 
               summary = FALSE)
fit2 <- fitted(m1b, 
               newdata=D2,
               #re_formula=NA, 
               allow_new_levels=TRUE, 
               #                   sample_new_levels = "old_levels", 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

sc_tmp_diff <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff > 0))

sc_tmp_diff2 <- sc_tmp_diff %>% 
  arrange(m) %>% 
  filter(p > .9)

cc <- cnl_dat %>% 
  mutate(country = countrycode(cowcode, origin="cown", destination = "country.name")) %>% 
  group_by(cowcode) %>% 
  summarise(country = first(country))

sc_tmp_diff2 <- left_join(sc_tmp_diff2 , cc)
sc_tmp_diff2 %>% mutate(obs = (1:n())/n()) %>% filter(cowcode %in% c(750, 70))  

variable <- "cspi_w"
D1 <- D2 <- cnl_ctrl_dat %>% 
  mutate(tv = as.numeric(tv))
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- -.5*s
D2[[variable]] <- .5*s

fit1 <- fitted(m1bc, 
               newdata=D1,
               #re_formula=NA, 
               allow_new_levels=TRUE, 
               sample_new_levels = "gaussian", 
               summary = FALSE)
fit2 <- fitted(m1bc, 
               newdata=D2,
               #re_formula=NA, 
               allow_new_levels=TRUE, 
               sample_new_levels = "gaussian", 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

sc_tmp_diffc <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff > 0))

sc_tmp_diff2c <- sc_tmp_diffc %>% 
  arrange(m) %>% 
  filter(p > .9)

cc <- cnl_dat %>% 
  mutate(country = countrycode(cowcode, origin="cown", destination = "country.name")) %>% 
  group_by(cowcode) %>% 
  summarise(country = first(country))

sc_tmp_diff2c <- left_join(sc_tmp_diff2c, cc)
sc_tmp_diff2c %>% mutate(obs = (1:n())/n()) %>% filter(cowcode %in% c(750, 70))  


variable <- "cspi_w"
D1 <- D2 <- cnl_dat %>% 
  select(cowcode, vp3_b, vp3_w, words, cspi_w, tv, spell) %>% 
  mutate(tv = as.numeric(tv)) %>% 
  filter(tv == 1)
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- -.5*s
D2[[variable]] <- .5*s

fit1 <- fitted(m2b, 
               newdata=D1,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m2b, 
               newdata=D2,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

sc_tmp_diffb <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff > 0))


sc_tmp_diff2b <- sc_tmp_diffb %>% 
  arrange(m) %>% 
  filter(p > .9)

variable <- "cspi_w"
D1 <- D2 <- cnl_ctrl_dat %>% 
  select(cowcode, vp3_b, vp3_w, words, cspi_w, tv, exconst_w, 
         gdppc_w, gdp_growth_w, war_w, EFindex_w, spell) %>% 
  mutate(tv = as.numeric(tv)) %>%
  filter(tv == 1)
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- -.5*s
D2[[variable]] <- .5*s

fit1 <- fitted(m2bc, 
               newdata=D1,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m2bc, 
               newdata=D2,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

sc_tmp_diffbc <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff > 0)) 

sc_tmp_diff2bc <- sc_tmp_diffbc %>% 
  arrange(m) %>% 
  filter(p > .9)

tmp_diff2_sc <- bind_rows(sc_tmp_diff2 %>% mutate(model = "Model 1"), 
                          sc_tmp_diff2b %>% mutate(model = "Model 2"), 
                          sc_tmp_diff2c %>% mutate(model = "Model 3")) %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name"))


ggplot(tmp_diff2_sc, aes(x=m, xmin = lwr, xmax=upr, y=reorder(name, m, mean))) + 
  geom_pointrange() + 
  geom_vline(xintercept=0, linetype=3) + 
  theme_bw() + 
  theme(panel.grid=element_blank(), 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="Change in Predicted Probability of\nAmendment Ratification", 
       y="", 
       caption="See Table A6 (models 1-3) for model results") + 
   facet_wrap(~model) 
  # ggtitle("Effect of Social Capital on Amendment Ratification")
ggsave("figures/sc_cnl_amendments_v2.png", height=12, width=8, units="in", dpi=300)


## Figure A4

variable <- "vp3_w"
D1 <- D2 <- cnl_dat %>% 
  select(cowcode, vp3_b, vp3_w, words, cspi_w, tv, spell) %>% 
  mutate(tv = as.numeric(tv))
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- D1[[variable]]-.5*s
D2[[variable]] <- D2[[variable]] + .5*s

fit1 <- fitted(m1b, 
               newdata=D1,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m1b, 
               newdata=D2,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

cr_tmp_diff <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff < 0)) %>% 
  arrange(m)

cr_tmp_diff2 <- cr_tmp_diff %>% 
  arrange(m) %>% 
  filter(p > .9)


cr_tmp_diff2 <- cr_tmp_diff2 %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name.en"))

variable <- "vp3_w"
D1 <- D2 <- cnl_ctrl_dat %>% 
  mutate(tv = as.numeric(tv))    
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- D1[[variable]]-.5*s
D2[[variable]] <- D2[[variable]] + .5*s

fit1 <- fitted(m1bc, 
               newdata=D1,
               #re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m1bc, 
               newdata=D2,
               # re_form = NA, 
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

cr_tmp_diffc <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff < 0)) %>% 
  arrange(m)

cr_tmp_diff2c <- cr_tmp_diffc %>% 
  arrange(m) %>% 
  filter(p > .9)


cr_tmp_diff2c <- cr_tmp_diff2c %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name.en"))

variable <- "vp3_w"
D1 <- D2 <- cnl_dat %>% 
  mutate(tv = as.numeric(tv)) %>% 
  filter(tv == 1)
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- D1[[variable]]-.5*s
D2[[variable]] <- D2[[variable]] + .5*s

fit1 <- fitted(m2b, 
               newdata=D1,
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m2b, 
               newdata=D2,
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

cr_tmp_diffb <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff < 0)) 

cr_tmp_diff2b <- cr_tmp_diffb %>% 
  arrange(m) %>% 
  filter(p > .9)

cr_tmp_diff2b <- cr_tmp_diff2b %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name.en"))

variable <- "vp3_w"
D1 <- D2 <- cnl_ctrl_dat %>% 
  mutate(tv = as.numeric(tv)) %>% 
  filter(tv == 1)
s <- sd(cnl_dat[[variable]], na.rm=TRUE)
D1[[variable]] <- D1[[variable]]-.5*s
D2[[variable]] <- D2[[variable]] + .5*s

fit1 <- fitted(m2bc, 
               newdata=D1,
               allow_new_levels=TRUE, 
               summary = FALSE)
fit2 <- fitted(m2bc, 
               newdata=D2,
               allow_new_levels=TRUE, 
               summary = FALSE)

fit1 <- as_tibble(t(fit1)) %>% 
  setNames(paste0("t", 1:nrow(fit1)))

fit2 <- as_tibble(t(fit2)) %>% 
  setNames(paste0("t", 1:nrow(fit2)))

tmp1 <- bind_cols(D1, fit1)
tmp2 <- bind_cols(D2, fit2)

tmp1 <- tmp1 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")
tmp2 <- tmp2 %>% 
  pivot_longer(starts_with("t"), 
               names_to = "iter", 
               values_to="vals")

cr_tmp_diffbc <- tmp1 %>% 
  ungroup %>% 
  mutate(diff = tmp2$vals - vals) %>% 
  group_by(cowcode) %>% 
  summarise(m = mean(diff), 
            lwr = unname(quantile(diff, .05)), 
            upr = unname(quantile(diff, .95)), 
            p = mean(diff < 0))

cr_tmp_diff2bc <- cr_tmp_diffbc %>% 
  arrange(m) %>% 
  filter(p > .9)

cr_tmp_diff2bc <- cr_tmp_diff2bc %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name.en"))

tmp_diff2 <- bind_rows(cr_tmp_diff %>% mutate(model = "Model 1"), 
                       cr_tmp_diffb %>% mutate(model = "Model 2"), 
                       cr_tmp_diffc %>% mutate(model = "Model 3"), 
                       cr_tmp_diffbc %>% mutate(model = "Model 4")) %>% 
  mutate(name = countrycode(cowcode, "cown", "country.name"))

m1 <- tmp_diff2 %>% 
  filter(model == "Model 1") %>% 
  rename(m1 = m) %>% 
  select(cowcode, m1)

tmp_diff2 <- left_join(tmp_diff2, m1) %>% 
  mutate(name = as.factor(name), 
         name = reorder(name, m1, mean), 
         lwr = case_when(lwr < -.06 ~ -.06, 
                         TRUE ~ lwr), 
         upr = case_when(upr > .06 ~ .01, 
                         TRUE ~ upr),
         credible = factor(as.numeric(sign(lwr) == sign(upr)), levels=c(0,1), labels=c("No", "Yes"))) %>% 
  filter(lwr != 0 & upr != 0 & name != "Iceland")

ggplot(tmp_diff2, aes(x=m, xmin = lwr, xmax=upr, y=name, color=credible)) + 
  geom_pointrange(show.legend = FALSE) + 
  geom_vline(xintercept=0, linetype=3) + 
  theme_bw() + 
  theme(panel.grid=element_blank(), 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="Change in Predicted Probability of\nAmendment Ratification", 
       y="", 
       caption="See Table A6 (models 1-4) for model results") + 
  facet_wrap(~model, ncol=2) + 
  scale_colour_manual(values = c("gray65", "black"))  
  #ggtitle("Effect of Constitutional Rigidity on Amendment Ratification")
ggsave("Figures/vp_cnl_amendments_v2b.png", height=10, width=8, units="in", dpi=1000)


## Figure 2
tmp_both <- cr_tmp_diff2c %>%
  select(-name) %>% 
  left_join(cc) %>% 
  mutate(eff = "Constitutional Rigidity") %>% 
  bind_rows(., sc_tmp_diff2c  %>% mutate(eff = "Social Capital")) %>%
  mutate(country = case_when(cowcode == 385 ~ "Norway", 
                             cowcode == 920 ~ "New Zealand", 
                             cowcode == 349 ~ "Slovenia",
                             grepl("Macedonia", country) ~ "Macedonia",
                             TRUE ~ country))


tmp_both <- tmp_both %>% 
  mutate(m2 = case_when(eff == "Constitutional Rigidity" ~ -m, 
                        TRUE ~ m)) %>% 
  group_by(country) %>% 
  mutate(oval = max(m2)) %>% 
  ungroup %>% 
  mutate(country = reorder(factor(country), oval, mean))




ggplot(tmp_both, aes(x=m, xmin = lwr, xmax=upr, y=country)) + 
  geom_pointrange() + 
  geom_vline(xintercept=0, linetype=3) + 
  theme_bw() + 
  theme(panel.grid=element_blank(), 
        plot.caption.position = "plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  labs(x="Change in Predicted Probability of\nAmendment Ratification", 
       y="", 
       caption="See Table A6 (model 1) for model results")  +
  facet_wrap(~eff)
  #ggtitle("Effect of Constitutional Rigidity on Amendment Ratification") 
ggsave("figures/fig2.png", height=12, width=7, units="in", dpi=1000)



## Table A7
cnl_refer_dat <- import("cnl_dat_refer.dta")
cnl_refer_ctrl_dat <- import("cnl_ctrl_refer_dat.dta")
m1br <- brm(amendment ~ vp3_b + words + cspi_w*referen + tv + 
              vp3_w:tv + (1 | cowcode) + (1|spell), 
            data=filter(cnl_refer_dat, lag_amend == 0), 
            family=bernoulli(),
            control = list(adapt_delta=.9), 
            cores = 4, 
            backend = "cmdstanr", 
            refresh=0, 
            silent=2, 
            seed=207, 
            save_pars = save_pars(all=TRUE))

m1bcr <- brm(amendment ~ vp3_b + words + cspi_w*referen + tv + 
               vp3_w:tv + exconst_w + gdppc_w + gdp_growth_w + EFindex_w + war_w + 
               (1 | cowcode) + (1|spell), 
             data=filter(cnl_refer_ctrl_dat, lag_amend == 0), 
             family=bernoulli(),
             control = list(adapt_delta=.9), 
             cores = 4, 
             backend = "cmdstanr", 
             refresh=0, 
             silent=2, 
             seed=207, 
             save_pars = save_pars(all=TRUE))

s3f <- summary(m1br)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,6,7,3,4,8,5,2), 
         p = unname(c(sigf(m1br)))) %>% 
  arrange(obs) %>% 
  select(-obs)

s3r_c <- summary(m1br)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m1br, "cowcode"))))

s3r_s <- summary(m1br)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m1br, "spell"))))

s3 <- bind_rows(s3f, s3r_c, s3r_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") 


s3$term <- ifelse(s3$term == "tv:vp3_w", "vp3_w", s3$term)


s3 <- s3%>% 
  select(-param)

s3$term[seq(2,20,by=2)] <- ""

s3 <-s3 %>% 
  mutate(term = gsub("vp3", "Const. Rigidity", term), 
         term = gsub("tv", "Time Varying", term), 
         term = gsub("cspi", "Soc. Capital", term), 
         term = gsub("judrev", "Jud. Review", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("referen", "Referendum", term), 
         term = gsub("sigma_", "Sigma: ", term), 
         term = gsub("_w", " (within)", term), 
         term = gsub("_b", " (between)", term))

#s3$vals2 <- ifelse(is.na(s3$vals), "", s3$vals)
s3 <- setNames(s3, c("Term", "Model 5"))

x0 <- cnl_dat %>% filter(lag_amend == 0)
x1 <- cnl_dat %>% filter(lag_amend == 0 & tv == 1)
dim(x0)
dim(x1)
length(unique(x0$cowcode))
length(unique(x1$cowcode))
aux <- data.frame(Term = c("N", "N Countries"), 
                  `Model 5` = c("2002", "80"))
names(aux) <- names(s3)
s3 <- bind_rows(s3, aux)



s3fc <- summary(m1bcr)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,6,7,3,4,8,9,10,11,12,13,5,2), 
         p = unname(c(sigf(m1bcr)))) %>% 
  arrange(obs) %>% 
  select(-obs)

s3rc_c <- summary(m1bcr)$random$cowcode %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "Country RE", term), 
         term = gsub("cspi", "Soc. Capital RE", term),
         p=unname(c(sigr(m1bcr, "cowcode"))))

s3rc_s <- summary(m1bcr)$random$spell %>% 
  select(c(1,3,4)) %>% as_tibble(rownames="term")  %>% 
  mutate(term = gsub("Intercept", "Spell RE", term),
         p=unname(c(sigr(m1bcr, "spell"))))

s3c <- bind_rows(s3fc, s3rc_c, s3rc_s) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="valsc") 

s3c$term <- ifelse(s3c$term == "tv:vp3_w", "vp3_w", s3c$term)


s3c <- s3c %>% 
  select(-param)

s3c$term[seq(2,30,by=2)] <- ""

s3c <-s3c %>% 
  mutate(term = gsub("vp3", "Const. Rigidity", term), 
         term = gsub("tv", "Time Varying", term), 
         term = gsub("cspi", "Soc. Capital", term), 
         term = gsub("judrev", "Jud. Review", term), 
         term = gsub("words", "log(# Words)", term), 
         term = gsub("EFindex", "Ethnic Frac.", term), 
         term = gsub("exconst", "Exec. Const.", term), 
         term = gsub("war", "Conflict", term), 
         term = gsub("referen", "Referendum", term), 
         term = gsub("gdppc", "log(GDP/capita)", term), 
         term = gsub("gdp_growth", "GDP Growth", term), 
         term = gsub("sigma_", "Sigma: ", term), 
         term = gsub("_w", " (within)", term), 
         term = gsub("_b", " (between)", term))

#s3c$vals2c <- ifelse(is.na(s3c$vals2c), "", s3c$vals2c)
s3c <- setNames(s3c, c("Term", "Model 6"))

z0 <- tmplbc %>% filter(lag_amend == 0)
z1 <- tmplbc %>% filter(lag_amend == 0 & tv == 1)
nrow(z0)
nrow(z1)
length(unique(z0$cowcode))
length(unique(z1$cowcode))

auxc <- data.frame(Term = c("N", "N Countries"), 
                   `Model 6` = c("1529", "78"))
names(auxc) <- names(s3c)
s3c <- bind_rows(s3c, auxc)

blnk_row <- c("", "")
s3new <- rbind(s3[1:16, ], blnk_row, blnk_row, blnk_row, blnk_row, blnk_row, 
               blnk_row, blnk_row, blnk_row, blnk_row, blnk_row, s3[17:22, ])
s3_all <- bind_cols(s3c, s3new %>% select(-1)) %>% select(1,3,2)

l1 <- loo(m1br, moment_match=TRUE)
l2 <- loo(m1bcr, moment_match=TRUE)

l_ests <- rbind(l1$estimates[3,], 
                l2$estimates[3,])


loos <- c("LOO IC (SE)", apply(l_ests, 1, function(x)sprintf("%.1f (%.1f)", x[1], x[2])))
names(loos) <- names(s3_all)
s3_all <- rbind(s3_all, loos)
r <- c("PRE", sprintf("%.2f", c(prefun(m1br), prefun(m1bcr))))
names(r) <- names(s3_all)
s3_all <- rbind(s3_all, r)
s3_all

